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A NUMERICAL METHOD FOR VARIATIONAL 
PROBLEMS WITH CONVEXITY CONSTRAINTS 

ADAM M. OBERMAN 



Abstract. We consider the problem of approximating the solution of 
variational problems subject to the constraint that the admissible func- 
tions must be convex. This problem is at the interface between convex 
analysis, convex optimization, variational problems, and partial differ- 
ential equation techniques. 

The approach is to approximate the (non-polyhedral) cone of convex 
functions by a polyhedral cone which can be represented by linear in- 
equalities. This approach leads to an optimization problem with linear 

I I constraints which can be computed efficiently, hundreds of times faster 

■^C than existing methods. 

-(— > 

rj Contents 

1. Introduction 2 

^ 1.1. Applications and related work 2 

^ 2. Counterexamples in approximating convex functions 4 

2.1. Failure of naive approaches 4 
Xi' 3. Polyhedral approximations to the cone of convex functions 5 
VjO 3.1. Directional convexity 5 
f~^ 4. Presentation of the variational problem 7 
^D 4.1. Example problems 8 
^^ 4.2. An analytic solution in one dimension 8 

• • 4.3. The Monopolist problem 10 

.J^ 4.4. A variant of the Monopolist problem 10 

^ 5. Approximation and implementation 11 

H 5.1. Overview of the discretization of the problem 11 

5.2. Convergence of convex approximations 11 

6. Numerical Discretization 11 

6.1. One dimensional constraint matrices 11 

6.2. One dimensional objective function gradient matrices 12 

6.3. Two dimensional objective function gradient matrices 12 

6.4. Two dimensional constraint function gradient matrices 13 

7. Accuracy of the convexity constraint 13 
7.1. Comparison of Wider Stencil Convexity Constraints 13 

8. Example Problems and Numerical Results 15 



Date: March 14, 2012. 



ADAM M. OBERMAN 

8.1. Projections d = 1 15 

8.2. Projections d = 2 15 

8.3. Numerical results for the Monopolist problem 15 

8.4. Numerical results for the variant of the Monopolist problem 15 

9. Reformulating polyhedral objective functions as linear constraints 15 

9.1. £°° and i^ projection 19 

9.2. Pointwise Constraints 20 

10. Performance improvement via Conic Programming 20 
References 20 



1. Introduction 

In this article we consider the problem of approximating the solution of 
variational problems subject to the constraint that the admissible functions 
must be convex. This is a numerical approximation problem at the inter- 
face between convex analysis, convex optimization, variational problems and 
Partial Differential Equation (PDE) techniques. 

In the theoretical setting, including a convexity constraint poses no ad- 
ditional difficulties, since the cone of convex functions is itself a convex set. 
However, in a computational setting, this problem has proven to be sur- 
prisingly challenging. First and foremost is the lack of a computationally 
tractable characterization of the cone of convex functions. Second, there 
are various mathematical difficulties which arise when working with approx- 
imations of convex functions. Convexity is not stable under perturbations: 
while strictly convex functions are still convex under small perturbations, 
nonstrictly convex functions are not. 

In this work we present a polyhedral approximation of the cone of con- 
vex functions. This approximation is computationally efficient in the sense 
that it can be represented by a number of linear inequalities which is small 
compared to the size of the problem. It can be used to build inner (strictly 
convex) or outer (slightly nonconvex) approximations to the cone of convex 
functions. 

Our methods are computationally efficient: the computational time to 
solve the problem with convexity constraints is comparably to a problem 
with simpler linear constraints. The increased efficiency is due to the re- 
duction in the number of constraints to enforce (approximate) convexity. 
The results are significantly (hundreds of times) faster, and generally more 
accurate compared to the method of [1], which is currently the most efficient 
method. 

1.1. Applications and related work. The earliest application of varia- 
tional problems with convexity constraints is Newton's problem of finding 
a body of minimal resistance [5]. There the convexity constraint arises as 
a natural assumption on the shape of the body; see [5] for a discussion and 
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also see [10]. A modern application is to mathematical economics [16, 12]. 
These problems can often be recast as the projection of a function (in the 
L^(r2) or H^(Q,) norm) onto the set on convex functions defined on the 
bounded domain 0,. Geometric applications include Alexandroff's problem 
and Cheeger's problem; see [11] for references. For a discussion of applica- 
tions to economics and history of this problem, we refer to [8]. 

There have been a few different numerical approaches to this problem, 
which rely on adapting PDE techniques to the problem at hand. Early 
work [4] using PDE-type methods, did not make assertions about the con- 
vexity (or approximate convexity) of the resulting solutions. Later work 
by [6] and [11] identified some of the difficulties in working with convex 
functions. These difficulties suggest that a straightforward adaptation of 
standard numerical methods is not possible. The introduction of a large 
number (superlinear in the number of variables) of global constraints was 
required in order to ensure discrete convexity. 

A recent work on the problem is [1] (see also [2]). In these works, approx- 
imate convexity is sought by enforcing positive definiteness of the discrete 
Hessian. However, as the authors of this work explain (see also §2), the 
fact that a discrete Hessian is non-negative definite does not ensure that 
the corresponding points can be interpolated to a convex function. The re- 
sulting optimization problem is a conic problem, which is generically more 
difficult to solve than one with linear constraints. However the number of 
constraints is less, on the order of the number of variables, in contrast to [6], 
in which the number of (linear) constraints grew superlinearly in the number 
of variables. 

In [9] approximations are performed in the space of bounded Hessians, 
and as a result, convexity may fail pointwise. 

A natural characterization of convex functions on scattered data can be 
found in Boyd and Vandenberghe [3, §6.5.5]: it uses the supporting hyper- 
plane condition. However it requires the introduction of new variables, and 
results in a very large number of constraint equations, one for each pair of 
data points. In [6], the variational problem from [16] was solved numerically, 
and convex envelopes were also computed. The supporting hyperplane con- 
dition for convexity is used to obtain discrete convexity constraints. This 
initially also leads to a quadratic number of constraint equations. By tak- 
ing advantage of the fact that points lie on a regular grid, the number of 
equations is reduced, but even after the reduction, the number of constraint 
equations is still superlinear in the number of grid points. 

More recently, [8] used the supporting hyperplane condition to derive 
convexity constraints. In this work the number of constraints is quadratic 
in the number of points used. In addition, the function is defined globally, 
essentially by the Legendre transform, so evaluation of function values and 
gradients is expensive. This idea of lifting hyperplanes to satisfy convexity 
conditions was also used in an early paper [15] on numerical solutions of the 
Monge- Ampere partial differential equation. 
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2. Counterexamples in approximating convex functions 

The difficulty in working with discrete approximations of convex func- 
tions is that the cone, C, of convex functions is a non-polyhedral cone with 
boundary. This leads to inequality constraints, which are more challenging 
to work with than equality constraints. In particular, equality constraints 
have been treated more frequently in the numerical PDE literature, e.g., the 
incompressibility constraint in the Stokes equation. 

The polyhedral convex functions are on the boundary of C: arbitrarily 
small perturbations of the these functions can be nonconvex. These func- 
tions, along with smooth but non-strictly convex functions, are the ones 
which cause difficulties (and these are the functions used to build coun- 
terexamples examples below). Strictly convex functions are generally easier 
to approximate, because for these functions, the convexity constraint is not 
active. So the challenge is to build approximations to convex functions which 
are robust on nonstrictly convex functions. 

2.1. Failure of naive approaches. Naive approaches to characterizing 
discrete convex functions can be shown to fail. The first (which naturally 
corresponds to a finite element approach) is to fix a triangulation and work 
with convex functions on the given triangulation. This approach fails be- 
cause it severely restricts the admissible convex functions to a subset of all 
convex functions. See [7] for details. 

The second approach (which is more naturally suited to finite differences) 
is to discretize the condition that the Hessian of a convex function is positive 
definite. This approach fails because the discrete Hessian (using finite dif- 
ferences) is a test which can give false positives (i.e., be positive definite on 
nonconvex functions) and false negatives (i.e., fail to be positive definition 
on convex functions). 

Example 1 (Testing convexity in coordinate directions is insufficient). It is 
well known that convexity in coordinate directions is not enough to ensure 
convexity. For example, the function u{x, y) = xy is linear (and thus convex) 
in each coordinate direction, but not jointly convex in x and y. In particular, 
the function restricted to the direction y = —x is concave. 

This last example shows that convexity must be enforced all directions. 
The approach of [6] was to enforce convexity in all directions available on 
the grid. This required a superlinear (approximately 0{N^'^)) number of 
constraints in the number of grid points A^. The resulting constraints were 
complicated to implement and computationally intractable. 

Example 2 (The discrete Hessian test fails). Two examples in [1] show that 
the discrete Hessian test can fail. The first is an example of a non-convex grid 
function, which nevertheless has a positive definite discrete Hessian. The 
second example is a piecewise linear function u{x, y) = max(3;, ax + by,y — c) 
which is convex but has a discrete Hessian which is not positive definite. 
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On the other hand, for smooth (C'^), strictly convex functions, it is true 
that (for smah enough h) the discrete Hessian is positive definite [1, Theorem 
4.2]. The theoretical drawback of this approach is that polyhedral functions 
are excluded from the class of admissible convex functions: instead smoothed 
and convexified versions are included. For example, in their computation of 
the Monopolist problem (see §4.3), the result is a smoothed approximation to 
the solution. Our approach has the advantage that piecewise linear functions 
on the grid are admissible and the numerically exact solution is obtained 
for the Monopolist example (see §8.4). The practical drawback of their 
approach is that it involves a conic optimization problem, to enforce the 
condition that the Hessian is positive definite. Our optimization problems 
has linear constraint, which is generally simpler than conic constraints [3]. 

Example 3 (Interpolations of convex functions may not be convex) . Given 
a convex function, the piecewise-linear triangulation of the function may 
not be convex. The following example comes from [6]. Consider u(x,y) = 
(x + y)^ on [0, 1]^. Then the piecewise linear approximation of u using the 
two triangles with common edge given by the line y = x is concave. 

It is also shown in [6] that, for fixed triangulations, enforcing convexity 
condition only on adjacent triangles results in an admissible set which is 
strictly smaller than C. 

3. Polyhedral approximations to the cone of convex functions 

In this section, we describe the polyhedral approximation to the cone 
of convex functions. In fact, we bracket the cone of convex functions by 
a family of interior cones, and a family of exterior cones, which we call 
uniformly convex and mildly nonconvex, respectively. In the limit as the 
directional resolution goes to zero, both of these cones should approach the 
cone of convex functions. 

The method we present uses finite difference discretizations of variational 
problems, along a local characterization of the convexity constraints. It 
has the flavor of PDE methods. Indeed, we use results and estimates first 
obtained in solving a PDE for the convex envelope [13, 14]. 

Our approach is to approximate the (non-polyhedral) cone of convex func- 
tions by a polyhedral cone which can be represented by linear inequalities. 
This approach leads to simple linear constraints: directional second deriv- 
ative constraints for the continuous problem, and sparse linear constraints 
for the discretized problem. 

We give a local characterization of the cone of convex functions. This 
results in a constrained variational problem where the number of constraints 
is linear in the number of variables. 

3.1. Directional convexity. The set, C, of convex functions is nonpoly- 
hedral, which is difficult to enforce numerically. We approximate this set 
by the polyhedral set of directionally convex functions, over a collection 
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of direction vectors. Define the directional resolution of a collection of d- 
dimensional unit vectors {vi},i = 1, . . . ,k as 



dtheta 



(1) 



d9 = max niincos 



-\w^v^ 



Note that dO measures the maximum angle between an arbitrary vector, and 
the direction vectors for the stencil. In the two-dimensional case, dO is half 
the maximum angle between direction vectors. In practise, accurate results 
are obtained with coarse directional resolution, see Figure 2. 

We can quantify the degree of non-convexity (or uniform convexity) al- 
lowed by the scheme with directional resolution dO. In the theorem below, we 
assume for simplicity, that the function is twice-differentiable. However, he 
results can be extended to continuous convex functions by mollification and 
application of Alexandroff 's theorem, as was done in [8]. Another approach 
is to use viscosity solutions, as was done in [13]. 



prop : approxConvex Proposition 1. Let u{x) he a twice- continuously differentiable function de- 
fined on M". Let V = {wj}f=i be a set of direction vectors, with directional 
resolution dO. Assume dO < 7r/4. // 

d^u 
(2) 312^0' i = l,...,k, 



Propl 



QuantConvexity (3) 



dvf 

then u is nearly convex, in the sense that 

Ai 



Prop2 



where Ai < 
addition 

(4) 

then 
I uconvex | (5) 

and u is convex. 



- > -tan^ide), 

An 

< An are the eigenvalues of the Hessian of u at x. Lf in 



d u] 2/ ,/,N ( d^u 

mm < --TT > > tan (at') max < -—^ 

vi&v dvf Vi€V dvf 



Ai > 



Remark 4. The conditions (2) are linear inequality constraints. By intro- 
ducing additional variables, the condition (4) can also be implemented as 
linear constraint. This is achieved by the linear inequality constraints 

d^u , 
0< 7< -^ < A 
dvf 

and 

7 > tan^{de)A. 

Strictness in the inequalities can easily be forced, for example by adding the 
term A — 7 to the objective function. 
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Proof. Let u{x) be the given twice-differentiable function, and suppose the 
minimum of Ai[n](x) occurs at x. Let wi be the eigenvector corresponding to 
Ai. Let 9 be the (positive) angle between wi and the nearest grid direction 
Vi. By (1), 9 < dO. Decompose Vi = cos9wi + sin^w, where w is a unit 
vector orthogonal to wi. Then compute 

(Pu 

(cos 6wi + sin 6w) D u (cos 9wi + sin 9w) 



dv: 



= cos 9wi D uwi + sin 9w D uw + 2 sin 9 cos 9wi D uw 

= cos 9\i + sin 9w D uw 

< cos 9\i + sin 6* An. 

In the computation above, we have used the fact that wi is an eigenvalue 
and w is orthogonal to wi to eliminate the cross term, and we have also used 
the estimate w D^uw < Xn- 

Now if (2) holds, the previous calculation yields 

— >-tai?{d9), 

An 

since 9 < d9, and (3) is established. 

Now, performing a similar calculation to the one above, but setting Wn 
to be the eigenvalue corresponding to An, and v^ to be the nearest grid 
direction to Wn, we obtain 

d^u 



, r, > COS 6'A„ + sin 9Xi. 
dvi 



Assuming (4) holds, we can combine the previous two inequalities with 
(4) to obtain 

cos^ 9Xi + sin^ 6* An > i&u^{d9) (cos^ 6* An + sin^ 9Xi) 

> tan^{9) (cos^ (9 An + sin^ 9Xi) 

using the fact that tan^ {d9) > tan^(0). Simplifying gives 

cos^ 9Xi > sin^ 9Xi 

Recall the assumption that 9 < d9 < vr/4, which means that cos{9) > sin(0). 
Thus, Ai > and (5) is established. D 

4. Presentation of the variational problem 
Let $7 be a closed bounded convex subset of M" , and 
C = {n : $7 — 7- M I n is convex in il.}. 
Consider the variational problem subject to the convexity constraint 

(6) inf I(u), where liu) = / f(x,u(x),Vu(x))dx, 

uecnic J^ 

and /C is a closed convex subset of a given space X = H^{Cl) ov X = L^(i7). 
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For the purposes of the numerical optimization, we require that for each 
fixed X, f is either a quadratic function of u and Vu, or / is convex and 
homogeneous of order one. In the first case, after expressing the approxima- 
tion of C by Hnear constraints, we arrive at a quadratic program (QP), and 
in the second case we arrive at hnear program (LP). Both are standard con- 
vex optimization problems, and can be solved by commercial or academic 
software packages [3]. 

If / is lower semicontinuous, coercive, and strictly convex on /C, the ex- 
istence and uniqueness of a minimizer directly follows from standard argu- 
ments. 

4.1. Example problems. 

Example 5. A typical example is given by 

If{u) = / -\Vu\^ + fudx 
Jn 2 

subject to n G C, and possibly with Dirichlet boundary conditions. 
Example 6. The Rochet-Chone example, [16], is given by 



Irc{u) = / -iul + Uy) + xux + yuy -udxdy 



1 

subject to n > 0, and u ^ C, with no explicit boundary conditions. 

Example 7. A standard problem is the projection (in some norm) onto the 
set of convex functions, which is given by the minimizer of the functional 

Iuo{u) = ll-U-Mollx, 

where X = H^{n) or X = L^{n), subject to n G C. 

Example 8. The problem of finding convex envelopes of a given function 
is also included (see [6]). Given uq G L'^{Q), take X = L^(ri) and 

luoiu) = \\u-uo\\L2(n) 

along with the additional constraint 

{u £ L (0) \ u < Uq a.e.}. 

The convex envelope can also be computed directly using a PDE [13, 14]. 

4.2. An analytic solution in one dimension. In this section we give 
an explicit non-trivial analytic solution in one dimension. The analytic 
solution is useful for verification of the numerical method. In this one- 
dimensional case, the solution is the convex envelope of the minimizer of the 
unconstrained problem, as was shown in [6, §4]. However, as was also shown 
in [6, §4], in higher dimensions this is no longer the case. 
Consider 

f^ 1 

min J[n]= / -u^{x) + f{x)u{x)dx 
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over functions in H^ with Dirichlet boundary conditions u{—l) = u(l) = 0. 
The convexity constraint becomes simply Uxx > 0. Here we restrict to the 
special case of a step function, 

j-c x<0 
I +c X > u 

for some c > 0. The minimizer will solve Uxxi^) = f{x) on some unknown 
set, and will solve Uxx = outside. Furthermore, it's easy to see that the 
set on which Uxx > will be of the form [a, 1], for some a > —1. Thus the 
candidate solutions will be convex, differentiable functions which are linear 
on [— l,a] and quadratic on [a, 1]. 

Description of the admissible set. By imposing the differentiability restric- 
tion at X = a, we can write candidates for the minimizer as a one-parameter 
family 

^ (7) u^i-) = h"^l! ^. forxE[-l,a] 

' ^ ^^ ^^ \c{x-b){x-l)/2 for2;G[a,l] 

where m = —c{a — 1)^/4, and b = (a^ + 2a — l)/2. Note that each u°'{x) is 
a local minimizer of the functional J. 

Evaluation of J[u]. Now compute J[ti''] as a function of a. The result, 
arrived using elementary integration, is: 

J[m"(x)] = --r^c^{a - lf{2,a^ + 10a - 1). 
48 

Now we can find a local minimum of J, by setting ■^J[u°'{x)\ = 0, to arrive 
at 

-c^{l-a){a^ + 2a-l) =0 
which gives the minimum at a = v2 — 1 . 

Verification of the Euler- Lagrange equation. The solution of the uncon- 
strained problem is 

(-f(x)x(x + l)/2 for2;<0 

111 Xi ^ \ 

[+f{x)x{x - l)/2 for x>0. 

The convex hull of the solution is obtained by putting 6 = into (7), which 
gives a^ + 2a — 1 = and agrees with the solution. 

The global minimizer, along with the solution of the unconstrained prob- 
lem is plotted in Figure 1. 
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Figure 1. Minimizer of the constrained problem (solid 
line) , and the minimizer of the unconstrained problem. 



figl 



monopolist 



monopolist var 



4.3. The Monopolist problem. The Monopolist problem of [16] is given 
by the variational problem 

min/[n] = / c|Vm| /2 — u(x) + Vu(x) • xdx 



for n = [0, 1]2 c M? subject to 



u>0. 



The exact solution is given in [16]. The quantity of interest in the econom- 
ics problem is the gradient map, which determines the optimal production 
choice of the monopolist as a function of the distribution of the population, 
which depends on two parameters. In this case, the gradient map is a com- 
bination of a constant, a linear functions, and nonlinear functions. These 
parts of the map correspond to a production of a fixed good for a large part 
of the population, and varying degrees of customization for the other parts 
of the population. 

4.4. A variant of the Monopolist problem. A variant of the Monopolist 
problem was studied in [12] and computed by [1] . It has the special feature of 
a linear objective function, which was easier to compute using the methods 
of [1]. The problem is given by 



min/[M] = / n(x) — Vm(x) • xdx 
uec ' ' Jn 



for n = [0, 1]2 c M^ subject to 



< Ux,Uy < 1, 



m(0) = 0. 
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The analytical solution of this problem is given by 

u{x,y) = Taa,x{0,x — a,y — a,x + y — b), a = -, 6=-(4 — v2), 

o o 

with optimal value 

^W = ^(6 + ^^)^.0549. 

5. Approximation and implementation 

5.1. Overvie^v of the discretization of the problem. Consider the 
functional 

^"(■w) = X] f^^^' u{xi),Vu{xi)) dx, 

i 

where Vn(a;j) is a finite difference approximation. The discretization of the 
gradients in the objective function is performed using standard centered 
finite differences. Significant performance differences resulted from different 
discretizations. 

The convexity constraints are implemented via a set of linear inequalities, 
at each point, which are directional second derivatives: 

u{xi) < '-^ '-, (i = 1, . . . , n), {j = l,...,k) 

where Vj is a collection of direction vectors. 

We also assume that the set /C can be represented by linear inequalities. 
The result is a quadratic (or linear) minimization, with linear constraints. 

5.2. Convergence of convex approximations. Proofs of convergence of 
approximations of variational problems with convexity constraints can be 
found in [6] and in [8]. In this article, we focus on the structure of the 
resulting discrete optimization problems produced by the approximation to 
the cone of convex functions. 

6. Numerical Discretization 

In order to have a well-posed (and accurate) convex variational problem, 
we need to carefully build the constraint matrices. We found that specific 
choices of these constraint and objective function matrices make a significant 
difference in the solution time for the variational problem. The choice of 
boundary conditions is also important. For simplicity, we present small 
examples of the matrices used to generate the computational results. 

6.1. One dimensional constraint matrices. In one dimension, convexity 
is enforced only at interior nodes. The convexity constraint operator D^x is 
given for n = 5 by: 

^ r 1 -2 1 

Z),^ = -2 1-2 10 

"- 1-2 1 
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6.2. One dimensional objective function gradient matrices. It is im- 



portant that the operator corresponding to 



in the objective function is 



a square symmetric matrix which is zero on constant functions. 

This is accomplished as follows. Start with the forward and backward 
finite difference operators, on the largest possible stencil. For example, in 
one dimension on a small grid. 



DH 



1 

h 



-1 1 

-1 











1 

-1 1 



The operator corresponding to u^ in the objective function is given by 



[Dt + (Dt 



{Dtr D 



Tn+ 



1 



1 
-1 







-1 

2 

-1 








-1 

2 
-1 





6.3. Two dimensional objective function gradient matrices. The two 

dimensional versions of the operator corresponding to |Viip can be con- 
structed as a square symmetric matrix in a similar way. 

The gradient squared objective function corresponds to a symmetric, pos- 
itive definite Laplacian operator. The right way to build this operator with 
finite difference matrices is to average the full forward and backward differ- 
ence operators in each grid direction. This ensures that the final operator 
is symmetric. (In fact corner values can be included by doing this with 
diagonal operators if this is desired) . 

So we build 



{D^c, + R 



yy> 



{Dt + D-+D+ + Dy], 



where the operators correspond to the one dimensional case above in each 
coordinate. Thus, in the n = 3 case, the symmetric positive Laplacian 
matrix is the 9x9 matrix 



-D 



xx+yy 



h^ 



2 


-1 





-1 

















1 


3 


-1 





-1 

















-1 


2 








-1 











1 








3 


-1 





-1 











-1 





-1 


4 


-1 





-1 











-1 





-1 


3 








-1 











-1 








2 


-1 

















-1 





-1 


3 


-1 

















-1 





-1 
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table :nbrs 



Figure 2. Neighbors for the width 1 and 2 scheme. Grid 
directions in the first quadrant for the width 1, 2, 3, 4 
schemes. 

Width (w) (Additional) direction vectors d9 = tan~^(l/?i;)/2 tan^^dO) 

A7 
.056 
.026 
.015 

Table 1 . Directions vectors in the first quadrant for schemes 

of width w. Values of d9 and tan^ d6. 



1 


(1,0) 


(1,1) 






.39 


2 


(2,1) 


(1,2) 






.23 


3 


(1,3) 


(3,1) 


(2,3) 


(3,2) 


.16 


4 


(1,4) 


(4,1) 


(3,4) 


(4,3) 


.12 



fig: schemes 



6.4. Two dimensional constraint function gradient matrices. If the 

gradient matrix appears in the constraints, we can afford to use higher ac- 
curacy centered differences away from the boundary, and we use forward or 
backward differences near the boundary. This which yields 



Dr 



1 
2h 



2 


2 











1 





1 











-1 





1 











-1 





1 











-2 


2 



7. Accuracy of the convexity constraint 

Proposition 1 gives a characterization of the approximation to the cone 
of convex functions. In this section we perform numerical tests to measure 
the accuracy of this approximation. 

7.1. Comparison of Wider Stencil Convexity Constraints. In this 
section we investigate the effect of enforcing convexity in more directions. 
For a fixed grid size of 21 x 21, we compared the effect of enforcing convexity 
in the four directions (horizontal, vertical, and the two diagonals) on the 9 
point nearest neighbor grid, and the additional 4 directions available on the 
wider grid. These grid directions are presented in see Figure 2 and Table 1. 
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For the purpose of comparison, the L^ projection was performed on various 
convex and non-convex test functions. 

For more examples, the worst case predictions of Proposition 1 were overly 
pessimistic. Indeed, for many functions we tested, the minimizer was the 
same up to numerical accuracy for the different methods. For the following 
functions, the minimizer was the same, for all width 1, 2, 3 schemes: 

sin(27rx), —x , (x — 3y) , |a; — 3y| . 

Notice that the schemes are invariant on convex functions, but may be 
(wrongly) invariant on slightly nonconvex functions. 

Next we present examples which were sensitive to the number of directions 
in which convexity was enforced. 

Results with a spiky nonconvex function. For other functions, the results 
with more directions differed strongly. For example, (always on [0, 1]^) 

exp(-30((a;-.5)2 + (y-.5)2)) 

gave results which converged as the stencil widened. The error, measured 
as the maximum value between the width 4 solution and the width 1,2,3 
respectively, was 0.02,0.01,0.003, for n = 21, and 0.02,0.01,0.005, for n = 
31. 

Results with a nonconvex quadratic. The next example uses the nonconvex 
quadratic function, (x^ — ay^)/2 rotated by an angle 9: 

(8) 5"''(^,y) = 

cos^ 6* + a sin^ 6* \ o // ^ n ^\ f acos'^ 9 + sin^ 9\ n 
X + ((1 — Qj cos 9 sm 9)xy + ( y , 



2 / V 2 

with a < 1, so that Ai = a, A„ = 1, and the eigenvectors are at an angle of 
9 from the coordinate axes. 

The worst case angles, and threshold values of a for detection of non- 
convexity are given by d9/2, a from Table 1. 

9 = 7r/8, a = -0.055, for width 1 

9 = 7r/12, a = -.0264, for width 2 

Taking a simple example, nonconvex quadratics, x^ — .5?/^ on [—1,1]^, 
rotated by 0, '7r/2, 7r/4 radians. The minimizer in this last case was the same 
for both schemes of width 1 and 2. 

Results were consistent with the predictions of Proposition 1. However 
values of = vr/S, a = — .5, yielded the same minimizers. 

For 9 = vr/S, a = —.1, the 4 direction scheme failed to see the non- 
convexity of the function and returned it as a minimizer, while the wider 
scheme gave something more convex. The results differed by .04. 

With 9 = vr/S, a = —.05 the solutions of width 1 and width 2 were the 
same. For width 3, the solution differed by .018. 
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8. Example Problems and Numerical Results 

8.1. Projections d = 1. In this section, we perform projections in dimen- 
sion d = 1. While the more interesting projections are in two dimensions, 
this example is useful for performance benchmarking and for visualization. 

The function to be projected is 

f{x) = sin(7ra;) 

on [-1,1]. 

Projections in L^, L^, L'^,H^, Hq, and H^ with gradient constraints of 
the function sin(7rx) are presented in Figure 3. When needed, the constraint 
u{0) = was also included. 

8.2. Projections d = 2. In this section, we compute for comparison pur- 
poses the projections from [6] and [1]. These are projections is various norms 
with convexity constraints. The function to be projected is given by 

spiky] (9) /(x,y) = -(4 + 5:Ey2)exp(-30((x-.5)2 + (y-.5)2)). 

Plots of the original function and the L^,L^,L°°, //^, and Hq projections 
are presented in Figure 4. 

8.3. Numerical results for the Monopolist problem. The example of 
the Monopolist problem of §4.3 was computed. The solution, along with a 
histogram of the gradient map is given in Figure 6. The numerical solution 
captures the fact that the density of the map is concentrated at point, along 
a line, and then spread out over a quadrant. 
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8.4. Numerical results for the variant of the Monopolist problem. 

The example of the Monopolist problem of §4.4 is well suited to our methods, 
since our representation of convex functions includes the piecewise linear 
functions aligned with the coordinate and diagonal directions, which is the 
case for this example. We obtain the exact solution up to interpolation 
error, and the computational time is which is hundreds of times less than 
in [1]. 

The contour lines of the numerical solutions are presented in Figure 6. At 
this resolution, the numerical solution is indistinguishable from the exact 
solution. (This is in contrast to the solution in [1], where the contour lines 
were curved). The computation was performed using piecewise constant 
quadrature for the integral, and the trapezoidal rule. In both cases only the 
most narrow (width 1) stencil was used to enforce convexity. In both cases, 
the piecewise linear solution was captured to within interpolation error. 

9. Reformulating polyhedral objective functions as linear 

constraints 

In this section we review a technique for reformulating various convex 
optimization problems. This is needed to translate the optimization problem 
into standard forms (Quadratic and Linear Programs). This reformulation is 
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-0.5 0.5 

(a) L^ projection 




-0.5 0.5 

(c) L°° projection 




-1 -0.5 0.5 

(e) H^ projection 




-0.5 0.5 

(b) L^ projection 




-1 -0.5 0.5 

(d) Hq projection 




-1 -0.5 0.5 1 

(f) H^ with gradient constraints 



Figure 3. Projection in various norms of sin(7rx) onto the 
cone of convex functions. 
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necessary for the (faster) MOSEK implementation. On the other hand, the 
CVX implementation recognizes convex optimization problems, and does 
not require reformulation. 
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(a) The function 






(b) L^ projection 






(c) L^ projection 






(d) I/°° projection 






(e) H^ projection 








(f) Hq projection 



Figure 4. Projection in various norms of the function given 
by (9) onto the cone of convex functions. 
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We also observed that different formulations of the optimization problems 
can lead to varying computational performance. 
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Figure 5. Solution of the Monopolist Problem: contour plot 
of the solution, and histogram of the gradient. 
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Figure 6. Contour plot of the numerical solution of the Mo- 
nopolist Problem, n = 21. 
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Method \ n 
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16 


32 


64 
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Zeroth Order Quadrature 
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.010 
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Trapezoidal Rule Quadrature 


.05 


.005 
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.0008 


Method from 1 
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.033 


.017 


X 


X 



Table 2. Monopolist Problem: error {L°°) for the methods 
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We start with a method for writing L°^ and L^ norms (of either gradients 
or the variables) as a linear program (i.e. linear constraints and linear 
objectives). This is a standard technique see for example, see [3] and [1]. 

9.1. ^°° and i^ projection. The reformulation of the objective function 
maxfc \uk — gk\ is accomplished by setting the objective function to t and 
enforcing, t = max \u — g\, using the linear constraints 

dk — t < Uk < Qk + t-, for each k. 

Similarly, the reformulation of the objective function ^^ \uk — gk\ is ac- 
complished by setting the objective function to ^j. t^, and enforcing \uk — gk\ 
tk using the linear constraints 

gk - tk <Uk<gk + tk, for each k. 
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9.2. Pointwise Constraints. Dirichlet boundary conditions 

u{x) = f{x), X G dD, 

can be implemented using equality constraints. 

For variational problems with one degree of freedom, (e.g., H^ projec- 
tion) one equality constraint should be included to force uniqueness of the 
minimizer, for example by setting 

u(0) = 0. 

The latter constraint is not necessary for the problem to be well-posed, but 
it improves the conditioning of the problem and can speed up the solver. 
For example, in the case of H^ projection in one dimension, with n = 2001, 
the solution time improved from 32.3s to 8.7s. 

10. Performance improvement via Conic Programming 

Conventional wisdom is that quadratic programming is generally faster 
than conic programming, and this is certainly the case for the L^ projections. 

However, using QP, the H^ projections were much slower than the L^ pro- 
jections. On the other hand, [1] used a conic reformulation and found these 
solutions times of the H^ and L^ projections to be comparable (although in 
both cases slower than ours). 

The conic reformulation replaces the objective 2^(0^ Dx)u with the new 
variable t and the constraint t < (Dxu)^. This formulation takes advantage 
of a matrix factorization 

rp 
^XX ^x ^X' 

Our implementation of the QP did not take advantage of the matrix factor- 
ization, which may explain the performance difference. 

We reformulated the quadratic objective as a conic constraint. Computed 
this way, the H^ projection was no more costly than the L^ projection. The 
results are shown in Tables 4-6. 
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CP time 



QP time 
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1.4 



1001 2001 
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1.5 I 1.7 
Table 4. Run time for d = 1, L^ projection, using Qua- 
dratic Program and Conic Program. 
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Table 5. Run time for d = 1, H^ projection, using Qua- 
dratic Program and Conic Program. 



CP time 



QP time 



1.6 



1.4 



1.7 



1.5 



16 20 32 45 64 90 128 



1.9 



2.3 



2.3 



6.4 



5.5 13 



117 1090 



42 



101 



251 



'able 6. Run time for d = 2, H^ projection, using Qua- 
dratic Program and Conic Program. 
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